package piechowicz.szymon.fft;

// Java port of http://code.google.com/p/fft-cpp/
public class FFT
{
    
    private final int n, lgN;
    private final boolean inverse;
    private final Complex[] omega;
    private final Complex[] result;
    private final Complex current_omega, t, u;
    
    public FFT(final int n, final boolean inverse)
    {
        this.n = n;
        this.inverse = inverse;
        this.result = new Complex[n];
        for (int i = 0; i < n; ++i)
            result[i] = new Complex();
        current_omega = new Complex();
        t = new Complex();
        u = new Complex();
        
        int lgN = 0;
        for (int i = n; i > 1; i >>= 1)
        {
            ++lgN;
            assert((i & 1) == 0);
        }
        this.lgN = lgN;
        omega = new Complex[lgN];
        int m = 1;
        for (int s = 0; s < lgN; ++s)
        {
            m <<= 1;
            if (inverse)
                omega[s] = Complex.exp(new Complex(0, 2.0 * Math.PI / m));
            else
                omega[s] = Complex.exp(new Complex(0, -2.0 * Math.PI / m));
        }
    }

    public Complex[] transform(final Complex[] buf)
    {
        bitReverseCopy(buf, result);
        int m = 1, half_m;
        for (int s = 0; s < lgN; ++s)
        {
            half_m = m;
            m <<= 1;
            for (int k = 0; k < n; k += m)
            {
                current_omega.set(1, 0);
                for (int j = 0; j < half_m; ++j)
                {
                    t.set(current_omega);
                    t.multiply(result[k + j + half_m]);
                    u.set(result[k + j]);
                    result[k + j].set(u.re + t.re, u.im + t.im);
                    result[k + j + half_m].set(u.re - t.re, u.im - t.im);
                    current_omega.multiply(omega[s]);
                }
            }
        }
        if (inverse == false)
            for (int i = 0; i < n; ++i)
                result[i].divide(n);
        return result;
    }

    public static double getIntensity(final Complex c)
    {
        return c.abs();
    }

    public static double getPhase(final Complex c)
    {
        return c.arg();
    }
    
    private void bitReverseCopy(final Complex[] src, final Complex[] dest)
    {
        for (int i = 0; i < n; ++i)
        {
            int index = i, rev = 0;
            for (int j = 0; j < lgN; ++j)
            {
                rev = (rev << 1) | (index & 1);
                index >>= 1;
            }
            dest[rev].set(src[i]);
        }
    }
    
}
